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ABSTRACT 

We study the dynamical stability and fates of hierarchical (in semi-major axis) two-planet systems 
with arbitrary eccentricities and mutual inclinations. We run a large number of long-term numeri¬ 
cal integrations and use the Support Vector Machine algorithm to search for an empirical boundary 
that best separates stable systems from systems experiencing either ejections or collisions with the 
star. We propose the following new criterion for dynamical stability: a ou t(l — e ou t)/ [ai n (l + ei n )] > 

2.4 [max(/ii n ,/Tout)] 1 ^ 3 (aout/din) 1 ^ 2 + 1-15, which should be applicable to planet-star mass ratios 
/x; n ,/r ou t = 10“ 4 — 1(T 2 , integration times up to 10 s orbits of the inner planet, and mutual incli¬ 
nations < 40°. Systems that do not satisfy this condition by a margin of > 0.5 are expected to be 
unstable, mostly leading to planet ejections if > p 0 ut, while slightly favoring collisions with the 
star for /ii n < /i out . We use our numerical integrations to test other stability criteria that have been 
proposed in the literature and show that our stability criterion performs significantly better for the 
range of system parameters that we have explored. 

Subject headings: planetary systems ■ planets and satellites: dynamical evolution and stability 


1. INTRODUCTION 

More than ~ 50 exoplanet systems discovered by ra¬ 
dial velocity (RV) surveys are known to harbor at least 
two planets, and many of them are in eccentric and well- 
separated orbits. The search for and characterization of 
these planets in either RV or transit surveys is gener¬ 
ally a time-consuming task, and having an-easy-to-use 
and accurate dynamical stability criterion is important 
to constrain either the existence of extra planets in the 
systems or the orbital configurations of already confirmed 
planets. 

Another motivation for searching for a stability cri¬ 
terion comes from theoretical studies in which plan¬ 
ets can have a variety of fates depending on the dy¬ 
namical stability of a planetary system. For exam¬ 
ple, instability can lead to the formation of f ree-floating 
planets through pl a net ej ections (e.g.. iSumi et alJf2011t 
IVeras fc Raymond l2012t i and planets reaching very 
nearly parabolic orbits can collide with or be tidally 
disrupted by the host star, becom ing a possible 
source of stellar metal pollu t ion ('e.g.. iSandauist et al.1 
l2002t iZuckerman et al.l 1200. 'll IVeras et al.l 1 201 3f~ Sim- 
ilarly, long-term stable and well-spaced planetary sys¬ 
tems can evolve secularly (with no orbital energy ex¬ 
change) to for m close-in planets by high-eccentricity 
migration (e.g.. iNaoz et aTTl2011b jWu fc Lit,hwickll20lii : 

iTevssandier et aJJ l2013t lPetrovichll2015H . A simple cri¬ 

terion to decide the fate of a planetary system based on 
its observed orbital configuration can help to constrain 
the most likely evolutionary path of different exoplanet 
systems without using expensive long-term V-body ex¬ 
periments. 

There is no analytic stability crit eria for arbitrary ec- 
centricities and/or inclinations (see lGcorgakarakosl l2008l 
for a review), while th e curren tly available (sem i-) ernpir- 
ical criteria (e.g., lHarringtonlll972l : lEggleton fe Kiseleva! 
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Il995t iMardling fc Aarsethll2001ll have generally not been 
tested in the planetary regime (in which one body 
contains almost all the mass of the system) or for 
the long timescales (up to ~ 10 6 — 10 8 orbits) dur¬ 
ing which two-planet systems can still becom e unstable 
(|Veras fc Mustilll 12013b [Petrovich et al.ll2014f) . 

In this study, we search for empirical criteria to de¬ 
cide whether a hierarchical two-planet system is likely to 
remain stable for long timescales or lead to either ejec¬ 
tions or collisions with the host star. We extend pre¬ 
vious numerical work (see 1 12.21) by considering a wider 
range of planetary systems, with planets in eccentric 
and/or mutually inclined orbits, and much longer evo¬ 
lution timescales. We also use, for the first time, the 
Support Vector Algorithm in the context of dynamical 
stability analysis, and fully detail our implementation. 

2. PREVIOUS WORK ON THE STABILITY OF 
TWO-PLANET SYSTEMS 

In this section, we briefly summarize the previous work 
on the stability of two-planet systems. We will use some 
of the stability criteria that have been proposed in the 
literature as benchmarks to compare to our results in 

ssu 

2.1. Stability of close two-planet systems with low 
eccentricities 

If the orbits of two planets are guaranteed to never 
cross, precluding collisions between planets or strong 
gravitational interactions, then they are said to be Hill 
stable. It has been shown that the conservation of angu¬ 
lar mome ntum and energy can constrain Hill stable tra - 
jectories dMarchal fe Bozisl 1 1 982t IMilani fe Nobilill 19881) . 

For low eccentricities (e < 0-1), the H ill stability crite¬ 
rion can be written as (e.g.. lGladmanlll99 Si) 

— >2.4( Atin + /W ) 1/3 + l, (1) 

Clin 

where a ou t (ain) and /ii n (/x Qu t) are the semi-major axis 
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and planet-to-star mass ratio of the inner (outer) planets, 
respectively. For reference, Equation m implies that two 
Jupiter-like planets orbiting a Sun-like star are Hill stable 
for aout ^ 1.30ain. Note that the Hill stability criterion 
does not discriminate mean-motion resonances. 

The Hill criterion gives no information about the long¬ 
term behavior of the system, and repeated interactions 
between planets in Hill stable orbits can still lead to ei¬ 
ther ejections and/or collisions with the star. The orbits 
that are protected against either ejections or collisions 
with the star are referred to as Lagrange stable. Also, the 
systems that fail the Hill criterion can still avoid having 
close approaches and be long-term stable (see discussion 

While there is no analytic criterion for Lagrange sta¬ 
bility, numerical studies show that the Lagrange sta- 
bilit y boundary lies close to the Hi l l stability bound - 
ary (jBarnes fc Greenberg! 120061 120071 : iDeck et, al.ll2012D . 
Based on the first-order mean-motion reson a nce overlap 
criteri on (lWisdomlll98d : iDuncan et al.lll989tj , IDeck et all 
(120131) studied the conditions that can yield chaotic be¬ 
havior in a two-planet system. These authors give the 
following criterion for the onset of chaos (which implies 
instability in their experiments) for two-planet systems 
in circular orbits: 

- < 1.46 (/li n + Mout) ^ + 1- (2) 

Gn 

For reference, from this criterion two Jupiter-like plan¬ 
ets are Lagrange unstable for a ou t 1$ l-25ai n . A 
numerical refinement of the chaotic zone boundary, 
wh ich sets the stability condi tion ab ove, is provided 
bv iMorrison fc Malhotral (|2015l b Also, iVeras fc Mustilll 
(|2013l) numerically studied the relation between the Hill 
and Lagrange stability boundaries for different eccentric¬ 
ities. 

Based on this previous work, hierarchical (or well¬ 
spaced, say aout > 2ai„) and coplanar two-planet systems 
with low eccentri cities are all expected to be long-term 
stable (e.g., IMarzarill2014D . Thus, the question of long¬ 
term stability in hierarchical two-planet systems should 
be focused on planets in eccentric orbits. 


2.2. Stability of hierarchical two-planet systems with 
arbitrary eccentricities 

There is no analytic criterion for the stability of hier¬ 
archical two-planet systems in eccentric orbittQ and most 
previous works rely on numerical experiments and/or 
heuristic approaches. We summarize some of the dy¬ 
namical stability criteria proposed for hierarchical two- 
planet systems. For consistency, we express each stability 
boundary in the form 


_ ^OUt(l ^OUt) ^ y 

ain(l T ^in) 


( 3 ) 


where ei n (e ou t) is the eccentricity of the inner (outer) 
planet and Y is a function of the initial orbital elements 


2 Extensions to the first-order resonance overlap criterion to 
eccentric orbits require taking into acco unt higher-order mean- 
motion resonances (e.g.. [Deck et a.0l2() 1 3 i. A calculation for non¬ 
zero (but still low) eccentricities considering only first-order reso- 
nances in the tes t -parti cle approximation has been carried out by 
IMustill fe Wvattl fgQT2 ). 


and masses. This choice is motivated by our results in 
1 14.1.11 where we find that the single parameter that best 
descri bes t h e stab ility bou ndary is r ap . 

(i) lEggleton fc Kiseleva] ( 1995 1 studied the stability of 
hierarchical triple systems with a wide range of masses 
and define a system to be n-stable if it preserves the 
initial ordering of the semi-major axes of the orbits and 
there are no escape orbits for 10 n orbits of the outer 
planet. The authors find an empirical condition for 2- 
stability using a set of N-body integrations, which in the 
planetary regime (p- ln , Mout <C 1) becomes: 


l-4/h n /: 


,1/3 Pout 


1 + Mout 


-1/3' 


( 4 ) 


The authors tested this criterion for p ln , p out > 0.01 and 
the following sets of orbital elements: prograde copla¬ 
nar orbits with either ej n £ [0,0.9] and e ou t = 0 or 
e 0 ut £ [0, 0.9] and ej n = 0, and circular orbits with mu¬ 
tual inclinati ons i m £ ] 0,18 0°]. 

(ii) iMardling fe Aarsethl (l200lh made an analogy be¬ 
tween the stability against escape in the three-body prob¬ 
lem and the stability against chaotic energy exchange 
in the binary-tides problem, and derived a semi-analytic 
stability criterion. We modify their criterion for coplanar 
and prograde orbits to express in the form of Equation 
© as 


' ap 


>4/ 


MA01 

crit 


= 2.1 


,(1 ^out) 
(1 + ein) 


(1 T Mout) 


1 + ^out 


(l-eont) 1 / 2 ] 


2/5 

(5) 


This criterion does not include a dependence on e; n and 
Pin, but the authors claim that it is valid for all eccentric¬ 
ities and masses of the inner body. Also, this criterion 
was proposed in the context of stellar clusters where, un¬ 
like our study, the mass ratios are not too different from 
unity. 

(iii) The Hill stability criterion by IMarchal fc Bozisl 
( 1981) can be written in the planetary regime as 
( Gladmanlll993D : 

r , p> V- sJ di^), (6, 

where S satisfies the implicit equation 


(Min + Mout/J 2 ) r / 2 \ 1 / 2 , /i \l /2 c 

3 Min (1 Gn) A Mout (f ^out) 0 

(Min + Mont) L 

MinMout 


1 2 


—1 — 3 4/3 - 


(Min + Mout) 


4/3 


= 0. 


( 7 ) 


The Hill stability cond ition has been extended by 
IVeras fc Armitagg (I2004D and iDonnisonl (120061 . I20li() 
to arbitrary mutual inclinations i m . Even though the 
Hill stability might not determine the long-term sta¬ 
bility of a two-plan et system (see i J2.ID , we will use 
it as a benchmark Barnes & Greenberg] I2006L 120071) . 
iKopparapu fc Barnea i 20101 ) numerically studied the re¬ 
lation between Hill and Lagrange stability and provided 
fitting expressions to determine the relation between 
these boundaries. Their results should be applicable 
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to planetary systems consisting of one terrestrial-mass 
planet and one much more massive planet with initial 
eccentricities less than 0.6. In this work, we focus on a 
complementary regime in which both the inner and the 
outer planets have masses much (at least ~ 30 times) 
larger than the Earth and therefo re do not attempt to 
compare our results with those bv lKonnarapu fe Barnesl 
(12010) _ 

(iv) iGiuppone et all (1201311 proposed a semi-empirical 
stability criterion for eccentric two-planet systems based 
on Wisdom’s criterion of firs t-order mean-motion reso¬ 
nance overlap (lWisdomlll980f l. The authors argue that 
the initial value of the relative longitudes of pericenter 
A w = Win + Gi n — (w ou t + G out ) can have a significant 
effect on the stability boundary, where Wi n (w 0 ut) and Gi n 
(G ou t) are the argument of pericenter and the longitude 
of the ascending node of the inner (outer) orbits, respec¬ 
tively. For the most conservative case with A vo = 180° 
the stability boundary is given by 


> 013 = 1 + 1.57 


2/7 . 2/7 ( «out 

Mi„ + Mout — 
\ W-in 


( 8 ) 


The authors also provide expressions for the case in which 
the ellipses are initially aligned (A® = 0), but an expres¬ 
sion for arbitrary values of Avj is not provided. 


3. NUMERICAL SIMULATIONS 

We run iV-body simulations of planetary systems con¬ 
sisting of a host star and two planets. 

We use the publicly available Bulirsch-Stoer (BS) in¬ 
tegration algorithm of MERCU RY6.2 with accuracy pa¬ 
rameter e = 1CT 12 ([ Chambers! ll999h . We justify the 
choice of this algorithm because we are mostly interested 
in the evolution of dynamically active systems, where 
planets experience close encounters, and BS handles close 
encounters better than the other integration algorithms 
in MERCURY6.2. We simulate the evolution for a max¬ 
imum time £ ma x given in units of the initial period of the 

inner planet P ln ^ = 2-77 ( Gm s /af n A , where m s is the 

mass of the central star and ai n ,i is the initial semi-major 
axis of the inner planet. The orbital elements are given 
in astrocentric coordinates and the typical conservation 
of energy and angular momentum are better than ~ 10~ 4 
and ~ 10~ 6 , respectively. We ignore the effects from gen¬ 
eral relativistic precession and tides in our calculations. 


3.1. Initial conditions and input parameters 

In Table 1, we summarize the input parameters, ini¬ 
tial conditions, and outcomes of the different simulations, 
which are all described in the following subsections. Our 
fiducial simulation is 2pl-fiducial. 

The ratios between the mass of the planet and that 
of the host for the inner and outer orbits, /q n and /i out , 
respectively, are chosen from a uniform distribution in log 
over the range [0.1 Mj/Mq, 10 Mj/Mq). The planets are 
treated as point masses, not allowing for planet-planet 
collisions. We note that the systems that would have 
planet collisions are expected to be unstable. 

In all simulations we start with a semi-major axis ra¬ 
tio that is uniformly distributed in a 0 ut/c3n G [3,10]. 
Thus, we generally exclude from our initial conditions the 
lowest-order mean-motion resonances p : p+q , withp = 1 
and q = {2,3,4} (a ou t/ain = {1.58,2.08, 2.51}), which 


can have a strong effect on the dynamics of the planetary 
system. Higher-order resonances have a weaker effect, as 
the strength of the resonant potential is proportional to 
e q . 

We draw the eccentricities of the inner and outer or¬ 
bits from a uniform distribution in [0, 0.9] and impose 
an upper limit to the eccentricity of the outer orbit 
e 0 ut < 1 — ain/<Jout to avoid a crossing of the initial or¬ 
bits. Note that all the orbits very close to this boundary 
become unstable and thus do not contribute significant 
information to the derived form of the stability bound¬ 
ary. 

For our fiducial simulation 2pl-fiducial we initialize the 
mutual inclinations i m between the inner and outer plan¬ 
etary orbits from a Rayleigh distribution with parame¬ 
ters <7j = 1°: the corresponding mean and median mutual 
inclinations are 1°.25 and 1°.17. In the simulations 2pl- 
inc-0 and 2pl-inc-20 we fix i m = 0 and i m = 20°, while in 
2pl-inc-rand we initialize i m from a uniform distribution 
in [0,80°]. 

The arguments of pericenter, the longitudes of ascend¬ 
ing node, and the mean anomalies are all drawn from a 
uniform distribution in [0,360°]. 

3.2. Dynamical outcomes 

We classify the different dynamical outcomes into the 
following categories. 

1. Two planets: two planets remain in the system for 
a time f max . Within this category, we distinguish 
the systems in which the initial semi-major axes 
of both planets have changed at a final time t max 
by less than 10%: \ai n j—a- m ,i\/ a in,i < 0.1 and 
|a 0 ut,/ — Oout^l /a ou t,i < 0.1. These systems have 
experienced only a small orbital energy exchange. 
In the complementary category at least one of the 
planets has changed its initial semi-major by 10% 
or more. 

2. Ejection: one planet is ejected from the system, 
which we define to happen when the planet reaches 
a distance from the central star > 100ai ni i. Such 
planets would almost certainly escape the system 
because at this distance the planet is either in an es¬ 
cape orbit (i.e., eccentricity > 1) or will most likely 
soon reach a escape orbit by energy perturbations 
from the inner planet. 

3. Stellar collision: one planet collides with the star. 
This is the only scale-dependent outcome because 
it depends on our definition of the ratio between 
the stellar radius and the initial semi-major axis 
R*/din,i- We use a fiducial conservative value for 
collisions of i?*/ai n ,i = 10 -4 , equivalent to placing 
the inner planet at = 46.5 AU for a solar-like 
sta r. W e study the effect of larger values of R+/a\ n j 
in THZl 

We treat the planets as point masses, not allowing for 
collisions between planets. However, for two Jupiter-size 
planets orbiting a Sun-size star the ratio ai n ,i /Rj in our 
fiducial simulation is ~ 10 5 , which is high enough that 
the rate of collisions between planets is expected to be 
very small compare d to the rate of ejecti ons or collisions 
with the star le.g.. iPetrovicli et al.ll2014ll . 
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TABLE 1 

Summary of simulated systems and outcomes 


Name 

&out / &in 

^in? ^out 

inc. 

[deg] 

Min? Mout 

[Mj/Mq] 

tmax 

-N S ys 

2 pi. with 

As <o.l 

2 pi. with 

As >o.i 

Ejection 

(uw > lo2 ) 

Stellar Coll. 

(uw = 10 ' 4 ) 

2pl-fiducial 

U(x; 3,10) 

U(x; 0,0.9) 

Ray(l) 

U(log x; -1,1) 

10 s 

3567 

2319 

8 

911 

329 

2pl-fid~4 

U(x; 3,10) 

U(x; 0,0.9) 

Ray(l) 

U(log x; -1,1) 

10 4 

3567 

2917 

390 

212 

48 

2pl-fid-5 

U(x; 3,10) 

U(x; 0,0.9) 

Ray(l) 

U(log x; -1,1) 

10 5 

3567 

2672 

258 

482 

155 

2pl-fid-6 

U(x; 3,10) 

U(x; 0,0.9) 

Ray(l) 

U(log x; -1,1) 

10 6 

3567 

2480 

80 

747 

260 

2pl-fid- 7 

U(x; 3,10) 

U(x; 0,0.9) 

Ray(l) 

U(log x; -1,1) 

10 7 

3567 

2369 

16 

874 

308 

2pl-inc-0 

U(x; 3,10) 

U(x; 0,0.9) 

1 

U(log x; -1,1) 

10 7 

2000 

1389 

31 

422 

158 

2pl-inc-20 

U(x; 3,10) 

U(x; 0,0.9) 

20 

U(log x; -1,1) 

10 7 

2000 

1422 

12 

402 

158 

2pl-inc-rand 

U(x; 3,10) 

U(x; 0,0.9) 

U(x; 0,80) 

U(log x; -1,1) 

10 7 

5000 

3319 

37 

1144 

500 


Note. Pi n j is the initial period of the inner planet. U{x]x minj^max) is the uniform distribution with x m [ n < x < # ma x and Ray(:r) 
is the Rayleigh distribution with parameter x. 



Fig. 1.— Fraction of systems with different dynamical outcomes 
as a function of the maximum integration time i ma x in units of the 
initial orbital period of the inner planet P[ n i . The dashed lines 
indicate the systems with two surviving planets for which at least 
one of the orbits has either changed its initial semi-major axis by 
> 10% at tmax (blue) or not (yellow). 

3.3. Results 

From Table 1, we observe that most systems (~ 65%) 
in our fiducial simulation 2pl-fiducial have two planets 
by the end of the simulation. Within this category over 
99% are in secularly stable orbits in the sense that the 
planets have experienced only small orbital energy vari¬ 
ations relative to their initial energies (the rms A ajai of 
the systems with |Aa|/aj < 0.1 is ~ 0.3%, where a* is 
the initial semi-major axis). 

The second most common outcome (~ 26%) is a sys¬ 
tem with one planet ejection, followed by a system with 
a stellar collision (~ 9%). 

The branching ratios into the different dynamical out¬ 
comes depend on various parameters, which we study 
next. 

3.3.1. Effect of the integration timescale t max 


In Figure |T| and Table 1 we show the evolution of 
the different outcomes in our fiducial simulation as a 
function of the integration timescale t max . In Table 
1 the simulation 2pl-ftd-x corresponds to 2pl-fiducial at 

t — inxp . 

^max 1 in,i* 

From FigureQ]and Table 1, we observe that the number 
of systems with two planets decreases as a function of 
time (or f max ) at the expense of increasing the number of 
ejections and collisions with the star, as expected. This 
decrease is most rapid for f max < 10 6 Pi n> i, after which 
time the fraction of systems with two planets (black line) 
shows a much slower decrease. For instance, from Table 
1 we see that the number of two-planet systems decreases 
by ~ 12.6% in going from 10 5 to 10 6 Pi nii , while it does 
so only by ~ 2.4% from 10 7 to 10 8 Pi ni i. 

From Figure [l] we observe that the fraction of systems 
with planets having significant variations in their semi¬ 
major axes (|Aa/di| > 0.1, blue dashed line) decreases 
rapidly from ~ 11% at t max = 10 4 Pi nj i to < 0.5% at 
fmax > 10' Pin, i- Thus, almost all the of systems with 
two planets that survive for more than 10 7 Pi nj i have ex¬ 
perienced small orbital energy variations relative to their 
initial values and might be regarded as secularly stable 
systems. 

In conclusion, our simulations show that there is lit¬ 
tle variation in the branching ratios of the dynamical 
outcomes after integrating the systems for longer than 
~ 10 7 Pi nj i. After this time the systems with two planets 
are essentially all in secularly stable orbits. 

3.3.2. Effect of varying R^/a i n ,i 

In Figure [H we show the fraction of systems in 2pl- 
fiducial with different outcomes as a function of the ratio 
between the stellar radius and the initial semi-major axis 
of the inner planet, -R*/ a in,i- 

We observe that the number of collisions with the star 
(green lines) increases with R*/ai U} i, as expected. For 
instance, the fraction of collisions using P*/ai n ,i = 10~ 4 
(ain,i ~ 50 AU for a solar-radius star) is ~ 9%, while this 
ratio increases to ~ 15% for R^/a^i = 10 -2 ((M n ,i ~ 0.5 
AU for a solar-radius star). 

From Figure^ we observe that the fraction of systems 
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R*/ ^in,i 

Fig. 2.— Fraction of systems in 2pl-fiducial with different dy¬ 
namical outcomes as a function of the the ratio between the stellar 
radius and the initial semi-major axis of the inner planets .R* /a- in i. 
The lower panel is a zoom-in of the fraction of systems experiencing 
stellar collisions. 

with two planets (black line) decreases only slightly (~ 
1%) by increasing i?*/ai n y from 10 -4 to 10 -2 , while the 
fraction of ejections decreases more significantly (~ 20%) 
for the same range of i?*/ai n y. 

In summary, varying R*/a- ln ^ mainly affects the ra¬ 
tio between ejection and collisions with the star, while 
the fraction of stable and unstable (ejections or colli¬ 
sions with the star) remains roughly constant. We shall 
use our conservative fiducial value of R*/a lni i = 10~ 4 
to determine the stability boundary in our subsequent 
analysis. 


3.4. Effect of the mutual inclination 

In Figure [3] we show the fraction of systems in 2pl- 
inc-rand with different outcomes for different bins of the 
initial mutual inclination i m . 

The figure shows that the fraction of ejections de¬ 
creases from ~ 0.23 for i m < 10° to ~ 0.16 for i m G 
[30°, 40°]. This decrease is marginally significant and 
might be related to the expected reduction in the time 
at which the planets experience close approaches when 
the orbits have higher mutual inclinations. For the same 
range of mutual inclination i m < 40° the fraction of stel¬ 
lar collisions does not show a clear trend. However, we 
observe a statistically significant decrease from ~ 0.09 at 
i m G [30°, 40°] to ~ 0.06 for i m G [20°, 30°]. 

As we start increasing the mutual inclinations from 
i m ~ 40° there is a clear and nearly monotonic increase 
in the rate of both ejections and collisions with the star. 
This behavior might be expected since larger values of 
i m > 40° can excite Kozai-Lidov eccentricity oscillations 
with large amplitudes, which can either decrease the peri- 
center distance to < R^/a^i producing stellar collisions, 
or simply increase the apocenter distance of the inner 
planet, promoting close encounters with the outer planet. 
As a consequence, the fraction of systems with two plan¬ 
ets decreases significantly from ~ 0.73 for * m G [30°, 40°] 
to ~ 0.55 for i m G [70°, 80°]. 

In conclusion, the main effect of increasing the mutual 
inclination from ~ 40° is the enhancement of the rate of 
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Fig. 3.— Fraction of systems in 2pl-inc-rand at f max = i 0 7 
as a function of the initial mutual inclination of the two planets. 
Upper panel: systems with two surviving planets. Middle panel: 
systems with one planet ejection. Lower panel: systems with one 
stellar collision. The error bars indicate the Poisson counting errors 
for each inclination bin. 

ejections and collisions with the star. As we increase the 
mutual inclinations from < 10° to ~ 30° — 40° there is a 
marginally significant decrease in the rate of ejections. 


4. SUPPORT VECTOR MACHINE (SVM) AND STABILITY 
BOUNDARY 

Our main goal is to find a stability boundary that best 
classifies the different outcomes and is simple enough 
(e.g., it has a small number of parameters) to allow for 
easy interpretation and use. We shall assess the per¬ 
formance of such a classification by its degree of “com¬ 
pleteness,” defined as the fraction of systems with true 
outcome A' that are correctly classified as A, or the ratio 
between the number of true positives and the number of 
true positiv es plus the number of false negatives (e.g., 
Ilvezid et al.ll2014l ). 

We use the Support Vector Machine (SVM) algorithm 
(e.g., IVap nikl Il996h to separate the « 2 P i = 1,2,..., JV 2 P i 
systems with two surviving planets from the i e j = 
1,2, ..,N e j systems with planet ejections and the 2 star = 
1, 2,.., TVstar systems with stellar collisions. 

We start by defining a set of parameters a. to de¬ 
fine a classification boundary. We assume that a 
is a simple function of the initial orbital elements 
{^out/diin; ^inj ^out; fm! and the masses {/q n , /r ou t}- 
For instance, we will define one set of parameters as 
a = [r ap , /q r ( ] with r ap defined in Equation ([3j and for 
each system i = 1,2,.., N syst we have a vector aq. 

We separate the classes by a hyperplane 

f{a) = /3 0 +f3-a t , 


(9) 
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where /3q and (3 are constants obtained using SVM. We 
define the separating function /(cc) such that f(a) > 0 
corresponds to systems with two planets (that is, stable 
systems), while f(a) < 0 could be either ejections or 
collisions with the star. We classify only two classes at 
the time: ejections from two surviving planets in H4.ll 
and stellar collisions from two surviving planets in 94.21 
For each system i = 1, 2, N syst we calculate /(aq) 
and we can formally define the completeness for each 
outcome as: 


/2pl = 

fe j = 

/star — 


/( a i2pi) > °| 

N 2p \ 

f( a U j) < 0[ 

iVej 

l/Kstar) < 0| 


( 10 ) 

( 11 ) 

( 12 ) 


where | - | is the cardinality of the set of systems and 
/ 2 pi,/ej,/star G [0,1]. Thus, a function that perfectly 
separates ejections (stellar collisions) from two surviving 
planets has / 2p i = / e j = 1 (/ 2 P i = /star = 1), while 
a conservative stability boundary would have f e j — 1 
(/star = 1) and significantly smaller / 2p i. 

We train the SVM classifier using the fitesvm pack¬ 
age from Matlab 2015a with standardized variables and 
a linear Kernel. For our fiducial simulation, we show the 
performance of each separation (i.e., / 2p i, / e j, /star in Ta¬ 
ble 2 and 3) using the same data as that in the training 
set. We have checked that the completenesses change by 
< 1% when a different set with ~1800 systems and sim¬ 
ilar initial conditions is used to test the performance of 
the classification. 

As discussed in 93.31 the number of systems with two 
surviving planets in our simulations is always larger than 
the number of systems with either ejections or collisions 
with the star. Thus, the SVM algorithm would naturally 
tend to classify the stable systems with a higher com¬ 
pleteness than ejections or collisions with star. Since we 
would like to have a boundary that separates each class 
with similar completeness (f 2p i ~ f e j and f 2p i ~ /star), 
we use a cost matrix in the SVM algorithm such that the 
cost of classifying a system into class X if its true class 
is Y is Nx/(Nx + Ny)- By doing so, we assign a higher 
penalty to misclassifying a class with a smaller number 
of systems. 

In practice, this arbitrary procedure works well for 
defining boundaries with similar completenesses and it 
mostly changes the offset of /(a) by a small amount 
relative to the case with equal misclassification costs. 
Finally, we note that by artificially promoting a bet¬ 
ter classification of systems with either ejections or col¬ 
lisions with the stars (smaller sample) at the expense of 
a poorer classification of stable systems, we expect to 
find a more conservative stability boundary in the sense 
that a smaller number of unstable systems are in stable 
regions. 


4.1. Separation of ejections and two surviving planets 

We start by separating the systems with two surviv¬ 
ing planets and from those with ejections because these 
classes dominate the branching ratios, and we leave the 
separation of stellar collisions and two planets for 94.21 


In Table 2 we show the completeness for different sepa¬ 
rating functions f(a) in different simulations (see Table 
1). We also include a set of previously proposed sta bility 
boundaries from Equations (j4|), ([5]), ©, and © in 92.21 
Similarly, in Figure 0 ] we show the distribution of the ra¬ 
tio between the number of systems with two surviving 
planets (solid black line) and ejections (red black line) 
and the total number of systems with either two planets 
or ejections for the stability boundaries above. The best 
criteria are those with values of / closest to unity. 

4.1.1. A sinqle parameter stability boundary: 

f(a) = /3 0 + di« 

We start by constructing a stability boundary that only 
depends on one parameter, using our fiducial simulation 
2pl-fiducial. We choose the parameter to depend on only 
the initial orbital elements and ignore the masses because 
without the orbital elements the masses cannot predict 
the fate of a planetary system. 

From Table 1, we observe that by setting a = r ap we 
obtain the function f(a) = r ap — 1.83 in our fiducial 
simulation. For this boundary we find completenesses 
of / 2p i ~ 0.86 and f e j ~ 0.87. Recall that by setting 
Aw = 180°, the parameter r ap becomes a measure of the 
minimum distance d m j n between two non-crossing copla- 
nar orbits and d m i n = ai n (l + ei n )(r ap — 1). Thus, the 
boundary can be rewritten as d m i n = 0.83-aj n (l + ei n ); in 
words, the boundary classifies a system as stable if ini¬ 
tially its orbits have a minimum distance that is at least 
~ 83% of the apocenter distance of the inner planet. 

One could calculate the initial minimum distance of 
the two ellipses for arbitrary values Aw and construct a 
stability boundary using a more precise measure of the 
closest approaches between the planets. However, this 
approach has a few shortcomings: 

1. the relative orientation of the orbits seems to have 
little effect of the performance of the stability 
boundary. In Table 2 we show that adding the ex¬ 
tra parameter cos(Ag7) does not increase the values 
of f 2p i and f ei . 

2. The resulting expression is too complicated to be 
of any practical use. 

Based on the arguments above, we will ignore the de¬ 
pendence on the initial relative apsidal angles Aw in our 
subsequent analysis. We are aware that in a case-by-case 
basis, the relative orientation can certainly make a dif¬ 
ference for the stability boundar y (see the discussion in 
94. 1.41 and iGiunpone et al.ll20 1 3l ). 

In summary, we argue that the best single-parameter 
stability boundary is r ap = 1.83 because of its simple 
functional form and the relatively high values of com¬ 
pletenesses it achieves, f 2p \ ~ 0.86 and / e j — 0.87. 

4.1.2. A two-parameter stability boundary: 

f{a) = /3o + /3i«i + /32«2 

Based on our findings above that the best single pa¬ 
rameter to describe the stability boundary is r ap , we fix 
oc\ = r ap and vary the functional form of a 2 to search for 
a two-parameter stability boundary that best separates 
stable systems from those with ejections in 2pl-fiducial. 

We start by including the dependence on the planet-to- 
star mass ratios /q n and p , ou t in f(a). Motivated by the 
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TABLE 2 

Summary of functions /(<*) found with SVM and from other works used to separate stable systems from systems with 

EJECTIONS 


Simulation 

/(«) 

/2pl 

/ej 

/star 

2pl-fiducial 



Tap — 1-83 


0.89 

0.92 

0.89 

2pl-fiducial 



r a P + 0.2 cos Aro — 1.83 

0.89 

0.92 

0.87 

2pl-fiducial 



rap - 4/df - 1.40 


0.91 

0.95 

0.80 

2pl-fiducial 



rap - 0.1^? - 1-83 


0.89 

0.92 

0.83 

2pl-fiducial 


T’ap 

— 2.7(/x in + /Xout) 1 ^ 3 — 

1.47 

0.90 

0.93 

0.89 

2pl-fiducial 



rap - 3.8/if/ 7 - 1.25 


0.91 

0.95 

0.79 

2pl-fiducial 


f"ap 

- 0.82/jV 3 (a out /a in ) - 

1.27 

0.93 

0.95 

0.69 

2pl-fiducial 


^ap 

2A h i/ 3 (“out/ain) 1/2 

- 1.15 

0.94 

0.96 

0.75 

2pl-fiducial 


T’ap 

3.2/if/ 3 (a ou t/ain) 1/3 

- 1.15 

0.93 

0.96 

0.74 

2pl-fiducial 


T’ap 

2.1/if/ 7 (a 0 ut/a in ) 1/2 

- 1.03 

0.94 

0.96 

0.75 

2pl-fid-5 


Tap — 

2 - 4 / i i 1 / 3 ( aout / a in) 1/2 

- 0.81 

0.99 

0.96 

0.57 

2pl-fid-6 


Tap 

2 ' 4 / i i 1 / 3 ( a °at/ain) 1/2 

- 1.01 

0.97 

0.96 

0.65 

2pl-fid- 7 


Tap — 

2-4/if/ 3 (clout /iJin ) 1 A 

- 1.09 

0.95 

0.96 

0.71 

2pl-inc-0 


Tap 

2-4/if/ 3 (liout/liin) 1//2 

- 1.15 

0.93 

0.96 

0.68 

2pl-inc-20 


Tap 

2-4/if/ 3 (clout/liin) 1 ^ 2 

- 1.15 

0.92 

0.95 

0.71 

2pl-inc-rand 


Tap 

2 - 4 ^i 1 / 3 ( a °at/ain) 1/2 

- 1.15 

0.92 

0.89 

0.57 

2pl-inc-rand (z m < 

20°) 

Tap 

2 - 4 /if/ 3 (^out/din) 1 ^ 2 

- 1.15 

0.92 

0.95 

0.71 

2pl-inc-rand (20° < im 

< 40°) 

Tap — 

2 ' 4 / i i 1 / 3 ( a oat/ain) 1/2 

- 1.15 

0.92 

0.95 

0.73 

2pl-inc-rand (40° < 2 m 

< 60°) 

Tap 

2'4/if/ 3 (“out/din) 1/2 

- 1.15 

0.92 

0.90 

0.56 

2pl-inc-rand (60° < i m 

< 80°) 

Tap 

2 - 4 /if/ 3 (“out/“in) 1//2 

- 1.15 

0.92 

0.80 

0.41 

2pl-fiducial 



r _ V EK95 
a P 1 crit 


0.92 

0.84 

0.87 

2pl-fiducial 



„ _ yMAOl 

a P 1 crit 


0.69 

0.98 

0.93 

2pl-fiducial 



r yHill 

7a P 7 crit 


0.80 

0.77 

0.81 

2pl-fid-4 



r yHill 

7a P 1 crit 


0.80 

0.84 

0.97 

2pl-fiducial 



„ yGMCl3 

ra P 7 crit 


0.63 

0.99 

0.99 


dependence of Hill’s stability criterion on the planet-to- 
star mass ratios, we test the performance of the following 

parameters: a 2 = Ah / 3 , Ah// an d (Min + Aw) 1/3 - From 
Table 2, we observe that the parameter that performs the 

best among these choices is a 2 = Ah/ 3 because it reaches 
the highest completenesses, / 2p i ~ 0.89 and f e j ~ 0.92 
compared to / 2p i ~ 0.86 and / e j ~ 0.87 for the one- 
parameter boundary. Also, we notice that incorporating 
the parameter a*/ 3 provides almost no improvement in 
the completeness relative to the single-parameter bound¬ 
ary f(a) = r ap — 1.83 and SVM assigns a small multi¬ 
plicative coefficient of /3 2 cs —0.1. Finally, a stability 
boundary using the parameter (A^in+A t out ) 1 ^ 3 marginally 
improves the performance of the boundary relative to the 
single-parameter expression, but it performs worse than 

simply using Ah/ 3 . 

We have tried different power laws of the form a if n and 
found for 9 = (1/2,1/3,2/7,1/4} the highest complete¬ 
nesses are reached with either 1/3 ~ 0.333 or 2/7 ~ 
0.286 (see Table 2). The performance is significantly 
(marginally) worse using 9 = 1/2 {6 = 1/4). Based on 
these results our method does not distinguish between 


the performance of a boundary with 9 = 1 /3 , which 
is expected from Hill stability (IGladmanl I1993T) and a 
bou ndary with 9 = 2/7, expec ted from resonanc e over - 
lap (lWisdomlll980H . No te that IMustill fe Wvattl (12012H 
extended the work by (|Wisdoml I1980H to planets with 
non-zero eccentricities and found a different power law 
with 9 = 1/5, which is not favored by our experiments 
relative to either 0 = 1/3 or 0 = 2/7. 

We experimented with various simple functional forms 

fl(ain,a ou t,ein,e out ) in a 2 = Ah/ 9 and found that by 
setting g = (ai n /a on t) 1 ' with v > 0 tends to increase the 
completeness relative to g = 1. In Table 1, we show the 
results of the stability boundaries for v = (1,1/2,1/3} 
and found the best results with v = 1/2. With this index 
the stability boundary is 

/ = r ap - 2.4AJi/ 3 (a ou t/ai„) 1/2 - 1.15, (13) 

2/7 

while for a a 2 = /i h ( g we find: 

/ = fap - 2.1^ n /7 (oout/«in) 1/2 - 1-03. (14) 

The boundaries in Equations m and m yield the 
highest completenesses in our heuristic search, / 2p i = 




















Petrovich 





Fig. 4.— Distribution of the ratio between the number of systems with two surviving planets (solid black line) and ejections (solid red 
line) and the total number of systems with either two surviving planets or ejections (i.e., all systems ignoring collisions with the star) as a 
function of different stability boundaries. The vertical dashed-dotted blue lines indicate the regions for which > 95% of the systems to the 
left (right) consist of ejections (two planets). Panel (a): single-parameter boundary / = r-,. p — 1.83 with r ap = a ou t(l — e 0 ut)/uin(l + e in) 

(see >1 1.1. 111 . Panel (b): two-parameter boundary / = r ap — 2.4/rV 3 (a ou t/a ; n ) 1 / 2 — 1.15 (see 44.1.211 . Panel (c): / = r ap — Wl 1 ^ 95 from 
lEggleton fe Kiseleval (119951 4 in Equation ([TJ. Panel (d): / = r ap — from IMardling fe Aarsethl 11200 11 1 in Equation (0- Panel (e): 

/ = 7a.p — frorn lGIadmanl (119931 4 in Equation (B). Panel (f): / = r ap — from lGiuppone et al.l (120134 in Equation ©. Note 

that the horizontal axes of panels (d), (e), and (f) are different from those in panels (a) through (c). 


0.94 and / e j = 0.95. We decided to stop here because the 
completenesses reached are close to unity. Also, note that 
we have experimented by adding an extra parameter ct 3 
with various functional forms and found only a marginal 
increase in the completenesses (~ 1 %), which are at the 
expense of a more complicated expression of /. 

In summary, we have found that the mass of the 
outer planet carries no information regarding the sta¬ 
bility against planetary ejections in our simulations and 
we only need to know the mass of the inner planet. We 
have found two stability boundaries with different power- 
laws for the inner planet-to-star /q n (Eqs. [13] and [14]'). 
which perform the best based both on the completeness 
they reach when separating ejections and stable systems 
and on their simplicity. 

4.1.3. Stability boundary and the maximum integration 
timescale. 

We have found that a simple function that separates 
stable from unstable systems in 2pl-fiducial is given 


by Equation m- Here we study the effect of the 
maximum integration time on this stability boundary. 
To do so we write the stability boundary as r ap = 

2.4/i 1 1 ^ 3 (a 0 ut/ain ) 1 ^ 2 +7 and see how 7 varies with f max . 

In Figure [5] we show our results for 7 as a function 
of t m ax- We show similar results in Table 2, labeled as 
2pl-fid-x (e.g., 7 = 1.01 for t max = 10 6 Pi ni i). From Table 
2 , we observe that the functional form above can well 
separate the stable systems from ejections (/ 2 P i,/ e j > 
0.95) just by changing 7 . 

From Figure [5l we observe that the required value of 7 
increases monotonically from ~ 0.6 for log (t m ax/Pin,i) = 
4 to 1.15 for log (t max /Pin,i) = 8. This increase is ex¬ 
pected because systems with larger r ap (fixing the masses 
and semi-major axes) should become unstable later. We 
also observe that the coefficient 7 increases more than ~ 
3 times more rapidly with time from log (f m ax/Pm,i) = 4 
to log (tmax/Pin.i) = 6 than for log (t max /Pin,i) > 6 (see 
linear fits). 
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Fig. 5.— Coefficient 7 from the stability boundary r ap = 

2 - 4 t‘ 1 1 n / 3 ( a out/am ) 1/2 + 7 as a function of the maximum integra¬ 
tion time tmax in units of the initial orbital period of the inner 
planet P[ n The red and blue dashed lines indicate a linear fit for 
log (tmax/ .Pin, i) S [4,6] and log (t m ax/Fin,i) e [6,8], respectively 

From Figure 0 we fit the evolution of 7 for 
log (i m ax/-Pin,j) > 6 and find that the stability bound¬ 
ary is given by 

r ap = 2 .yf ( V2 + 0.069log (+ 0 . 6 . 

\ ^in / \-*in ,i / 

(15) 

This stability boundary is valid for imax/.Pin,i = 10 6 — 
10 8 and given the slow variation of 7 with time, it sug¬ 
gests that only a small fraction of stable systems in 2pl- 
fiducial can become unstable in longer timescales. 

4.1.4. Misclassified systems 

Some of the stable systems (~ 6 %) are classi¬ 

fied as ejections because they initially satisfy r ap < 

2Aii 1 J^ t {a on t/a\ n ) 1 ^ 2 + 1.15. These systems tend to start 
with relatively aligned orbits: ~ 50% (~ 82%) start with 
cosAcu > 0.8 (cosAtu > 0). We checked that in some 
extreme cases, the system starts with r ap < 0 and avoids 
orbit crossing by starting with cos Am ~ 1 and engaging 
in a secular resonance. 

Similarly, some of the systems with ejections (~ 4%) 
are classified as stable because they initially satisfy r ap > 

2.4/r p ^(aout/ain ) 1//2 + 1-15. These systems tend to start 
with relatively misaligned orbits: ~ 40% (~ 76%) start 
with cosAru > 0.8 (cosAtu > 0 ). 

By adding the extra parameter cos A m to our sta¬ 
bility boundary we find ^”ap — ^.4/^ ou t (u 0 ut/^in) ^ T 
0 . 2 cosAtu + 1.1 and the completenesses increase only 
marginally from / 2p 1 ~ 0.94 and f e j ~ 0.96 to / 2p 1 — 0.9*4 
and f e j cs 0.97. 

In conclusion, some of the misclassification might be 
explained by the initial relative orientation of the el¬ 
lipses since the orbits that start with more aligned 
(misaligned) pericenters tend to be more stable (unsta- 
ble). These results a re consistent with the claims by 
iGiunnone et al.1 (|2013f l. However, the overall effect of 


Am only marginally improves the performance from our 
simpler stability boundary. 


TABLE 3 

Summary of functions f(a) found with SVM used to 

SEPARATE STABLE SYSTEMS FROM SYSTEMS WITH STELLAR 
COLLISIONS 


Simulation 
2pl-fiducial 
2pl-fiducial 
2pl-fiducial 
2pl-fiducial 
2pl-fiducial 
2pl-fiducial r 


T’ap 

?”ap 


ap 


/(“) 

/2pl 

fej 

/star 

Tap — 1-83 



0.89 

0.91 

0.89 

rap + 0.47 m 1 / 3 - 

1.93 


0.89 

0.91 

0.89 

r ap — 3.4M 0 / t — 

1.45 


0.91 

0.87 

0.91 

- 2.7(/i.in + /Xout)!/^ — 

1.58 

0.90 

0.93 

0.89 

- (“out/a 

in) 

1.18 

0.93 

0.83 

0.92 

^•^out ( a out/ttin 

)l/2 

- 1.15 

0.92 

0.83 

0.92 


4.2. Separation of stellar collisions and two surviving 

planets 

Following the same procedure as in § §4TT] and 14.1.21 
we search for a stability boundary that separates systems 
that experience stellar collisions from systems with two 
surviving planets in our fiducial simulation 2pl-fiducial. 

In Table 3, we show our results for some separating 
functions found using SVM and their corresponding com¬ 
pletenesses. Similarly, in Figure [ 6 ] we show the distribu¬ 
tion of the ratio between the number of systems with two 
planets (solid black line) and collisions (solid red line) 
and the total number of systems with either two planets 
or collisions for different stability boundaries, including 
those in Equations (QJ-([6]), and © from H2.21 

From Table 3, we observe that the single-parameter 
boundary using r ap is given b y f(a) = r ap — 1.83, which 
is identical to that found in N4.1.11 for separating ejec¬ 
tions from systems with two surviving planets. The 
completenesses using this function are / 2p i = 0.89 and 
/star = 0.89. 

As in H4.1.21 we include the dependence on the planet- 
to-star mass ratios /q n and /t ou t hr f{a), and test the 
performance of the following parameters: cr 2 = /t 1 / 3 , 
/t// 3 , and (/ri n T/itout) 1 ^ 3 - From Table 3, we observe that 

the boundary with a 2 = /i ir ] does not improve the per¬ 
formance relative to the single-parameter boundary (the 
completenesses are the same). 

By setting a 2 = (/q n + /Lout ) 1 / 3 we observe that there 
is a slight improvement since / 2p i increases from 0.89 
in the single-parameter boundary to 0.9 and the result¬ 
ing separating boundary is very similar to the one found 
for separating ejections from stable systems (see Table 

2). Finally, by setting a 2 = p, out we find that the per¬ 
formance improves more significantly and the complete¬ 
nesses are / 2p 1 = / sta r = 0.91. We tried other functional 
forms for the mass ratios and observed no improvements 
relative to a 2 = /q/ 3 . 

Similar to our procedure in H4.1.21 we add dependence 
on the semi-major axis ratio a ou t/«in- We find that the 

best separation is reached by setting a 2 = /+/ 3 (aout/adn) 
with completenesses of / 2p 1 = 0.93 and / sta r = 0.92 (see 
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Fig. 6 .— Distribution of the ratio between the number systems with two surviving planets (solid black line) and stellar collisions 
(solid green line) and the total number of systems with either two surviving planets or stellar collisions (i.e., all systems ignoring 
ejections) as a function of different stability boundaries. The vertical dashed-dotted blue lines indicate the regions for which > 95% 
of the systems to the left (right) consist of stellar collisions (two planets). Panel (a): single-parameter boundary / = r ap — 1.83 with 

r ap = a ou t(l - e ou t)/a in (l + e in ) (see i| 1.2| l. Panel (b): two-parameter boundary / = r ap - 2.4/r^(a out /a ill ) 1 / 2 - 1.15 (see j|4.2t . Panel 
(c): / = r ap — from lEggleton fc Kiseleva! 119951 1 in Equation (3J. Panel (d): / = r ap — from IMardling fe Aarsethl 120011 1 

in Equation (SJ. Panel (e): / = r ap — from lGladmanl 1199311 in Equation {BJ. Panel (f): / = r ap — from IGiuppone et al.l 

J2013ft in Equation ©. Note that the horizontal axes of panels (d), (e), and (f) are different from those in panels (a) through (c). 


Table 3). A slightly worse separation (/ 2p i = /star = 

0.92) is reached with a 2 = Mout (aout/ain) 1 ^ 2 , where the 
function reads 

/ = f ap - 2.4^(aout/ain) 1/2 - 1-15. (16) 

This function is identical to that in Equation 031, found 
to separate ejections from stable systems if we change 
Mout for /iin- This surprising result suggests that the 
long-term stability of the system against either ejections 
or collisions might only depend on max(/.q n , /x 0 ut)- Then, 
if an unstable system has > p, out , the most likely 
outcome is an ejection, while a collision with the host 
star is slightly more likely if fj, ln < /r ou t- Motivated 
by these findings we favor the separating function in 

Equation 1131) over r ap — 0.9/r p ^(a out /oi n ) — 1.18. Fi¬ 
nally, we have experimented with different exponents v 

in a = Mout( a out/flin) J/ and found no improvement rela¬ 
tive to v = 1/2. 


In summary, the mass of the outer (and not the inner) 
planet carries most of the information about the systems 
that collide with the star. We found that a good stability 
boundary for separating collisions from stable systems is 
given by Equation (HU), which is identical to the one 
found for separating ejections from stable systems when 
changing /z ou t for /q n . 

4.3. A criterion for stability against either ejections or 
stellar collisions 

In i J4.ll and M.2I we have found criteria for separating 
systems with ejections from systems with two surviving 
planets and systems with stellar collisions from systems 
with two surviving planets, respectively. In Figure [T] we 
show these two criteria by plotting the stability bound¬ 
ary against ejections in Equation dl3l) versus the stabil¬ 
ity boundary against collisions with the star in Equation 
(USD for the different outcomes in 2pl-fiducial. We note 
that these criteria differ only on the whether the mass of 
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Fig. 7.— Boundary found f or se parating stellar collisions from 
stable systems from Equation GSJ as a function of the boun dary 
found for separating ejections from stable systems Equation (|13[) 
for the different outcomes in 2pl-fiducial. The diagonal dot-dashed 
blue line indicate boundary in which both boundaries cross (i.e., 
Min = Mout)- 


the inner and the outer is used. 

We observe that most systems in regions in the positive 
quadrant (Equations [13] and [10] are both positive) are 
stable (very few red and green dots). This result suggest 
that we can combine Equations m and m to define 
a stability condition against either ejections or collisions 
with the star as: 

r ap > 2.4 [max{^ in ,/Xout}] 1/3 ( a out /a in ) 1/2 + 1.15.(17) 

This criterion yields completenesses of / 2 P i — 0.9 and 
/ej+star ^ 0.95. Since / 2p i < / e j+star this criterion is 
somewhat conservative in the sense that a smaller num¬ 
ber of unstable systems are in stable regions relative to 
the number of stable systems in unstable regions. 

Similarly, we observe from Figure [7] that most systems 
in the negative quadrant are unstable and by using the 
minimum instead the maximum in Equation (11711 we get 
/ 2 P i ~ 0.81 and / e j+ s tar — 0.97. Moreover, within the un¬ 
stable systems we observe that almost all of the collisions 
with the star (~ 95%) are in regions where /q n < /x ou tj 
while most (~ 72%) of the systems with ejections have 
/x> /Xout- Since we have an overall higher rate of ejec¬ 
tions than stellar collisions, we find that both rates are 
comparable in regions where /ii n < p, ou t: — 45% and 
~ 55% of the unstable systems undergo ejections and 
collisions with the star, respectively. 

In summary, we combine our previous results in Equa¬ 
tion (ffZD to propose a stability boundary against either 
ejections or collisions with the star. Systems that are 
unstable and have /q n > /x out will most likely undergo a 
planet ejection, while the systems that have /.q n < fi ou t 
will have a similar rate of ejections and stellar collisions, 
with the latter being slightly higher. 

4.4. Effect of mutual inclinations on the stability 
boundary 


£ 
<D 
-t-3 
CO 



^ap 2.4 [max (/iin? Mout)] / (^out/^in) ^ 1.15 



Fig. 8 .— Fraction of systems with two planets (solid black line), 
either ejections or collisions with the star (solid blue line), ejections 
(dashed red line), and collisions with the star (dashed green line) 
in the simula tion 2pl-inc-rand as function of the stability boundary 
in Equation G3- Panel a: systems with mutual inclinations Zm 
40°. Panel b: systems with mutual inclinations z m < 40°. The 
vertical dashed gray lines indicate the regions for which > 95% of 
the systems to the left (right) consist of either ejections or collisions 
with the star (two planets). 


From Table 2 we observe that the stability bound- 
ary against ejections in Equation (fl3l) found using 2pl- 
fiducial performs relatively well in 2pl-inc-rand (/ 2p i — 
0.92 and / 2p i — 0.89), which has a random distribution 
of the mutual inclination in i m [0, 80°]. 

By taking different bins of i rn in 2pl-inc-rand we find 
that the performance of the boundary in Equation m 
is the same (/ 2p i — 0.92 and / 2p i — 0.95) for the systems 
starting with i m < 20° and i m G [20°, 40°]. Similarly, the 
performance of this boundary in the coplanar case (i m = 
0 ) 2pl-inc-0 and in the simulation 2pl-inc-20 with i m = 
20° is almost the same. Thus, our stability boundary 
against ejections performs well for mutual inclinations 
im < 40°. 

As we increase the initial mutual inclinations in 2pl- 
inc-rand we find that f e j drops from ~ 0.95 for i m < 40° 
to 0.9 and 0.8 for i m G [40°,60°] and i m G [60°,80°], re¬ 
spectively. On the contrary, / 2p i remains equal to ~ 0.92 
for all bins in mutual inclinations. As discussed in 43.41 
this behavior might be expected since larger values of 
i m > 40° can excite Kozai-Lidov eccentricity oscillations, 
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which can promote close encounters with the outer planet 
and produce ejections in regions that would be long-term 
stable for i m < 40°. 

In Figure 01 we show the fraction of systems with 
different outcomes in our simulation 2pl-inc-rand for 
i m < 40° in panel (a) and * m > 40° in panel (b) as 
a function of the stability criterion against either ejec¬ 
tions or stellar collisions in Equation m- From panel 

(a) we observe that there in only a small fraction of 
systems with either ejections or collisions with the star 
for r ap > 2.4 [max{^ in ,/i out }] 1/3 (a 0 ut/ain) 1/2 + 1-15 and 
fe j — /star — /star+ej - 0.96 for i m < 40°. From panel 

(b) we observe that this fraction of unstable systems in 
stable regions increases for z m < 40° and the complete¬ 
nesses decrease significantly: f e j ~ 0.87, / sta r — 0.67, 
and /star+ej — 0.8. Since / s t ar is significantly lower than 
f e j we conclude that the performance of our stability cri¬ 
terion worsens mostly at the expense of having collisions 
with the star in regions classified as stable. This effect is 
observed in Figure [ 8 ] as an increase in the tail with posi¬ 
tive value of Equation (da of the distribution of collisions 
with the star (green dashed line) in panel (b) relative to 
panel (a). 

In summary, our stability criterion in Equation da 
performs well for mutual inclinations i m < 40°. For 
higher mutual inclinations the Kozai-Lidov mechanism 
produces a significant fraction of unstable systems in re¬ 
gions classified as stable and our criterion becomes a poor 
predictor for long-term stability. Our resul ts seem con - 
sistent with a previous study bv iGeorgakarakosl ( 2013 1. 
which shows that the mutual inclination has very little 
effect on the stability boundary for i m £ [0,40°]. 

5. DISCUSSION 

The main results of this paper are a set of new stabil¬ 
ity boundaries for separating systems that become un¬ 
stable against ejections and collisions with the star from 
systems that retain their two planets with small orbital 
energy variations (Equations [13] and [TB] )■ In particu¬ 
lar, we propose that hierarchical two-planet systems are 
long-term stable if they satisfy the condition in Equation 

GZD- 

Additionally, we find that our stability boundary: 

1 . performs significantly better than other previously 
proposed criteria (see completenesses in Table 2 
and 3, and Figures 0]and ©; 

2. performs well for all mutual inclinations * m < 40°; 

3. and changes slowly with the maximum integration 
timescale as oc 0.071og(f max /P in ) for f max /P in = 
10 6 — 10 s , while it does so ~ 3 times more rapidly 
for tmax/Pin = 10 4 —10 6 (see Figure[5]and Equation 

my 

The fate of the unstable systems depends mostly on the 
planetary masses. Most systems with /q n > /i out lead to 
ejections, while for fi m < p out there is a slightly higher 
number of collisions with the star than ejections. 

In what follows we discuss some of the consequences 
of our findings in the context of other works and the 
observations. 


5.1. Performance of other stability criteria 

In Table 2 we show the completenesses of the differ¬ 
ent stability criteria discussed in 42.21 that were reached 
in our fiducial simulation 2pl-fiducial. Similarly, panels 
(b) to (f) in Figures 0] and [ 6 ] show the fraction of out¬ 
comes for each stability boundary. In these figures the 
performance is best when the outcomes have the sharpest 
transition from 0 to 1 (a perfect separation leads to two 
step-functions). 

First, our simulations show that the Hill stability cri¬ 
terion (Equation [5] for coplanar systems) is a poor 
indicator of the dynamical stability of two-planet sys¬ 
tems because it achieves relatively low completenesses 
(/2pi ~ fe j ~ /star ~ 0.8). For comparison, even the 
single-parameter boundary r ap = 1.83 performs signifi¬ 
cantly better (/ 2 P i ~ f e j ~ /star ~ 0.9). Moreover, an im¬ 
portant fraction of the systems that are classified as Hill 
unstable are actually long-term stable (see the solid black 
lines in panel (e) of Figures 0] and [ 6 ] with r ap < F™ 1 ). Q 

Second, we observe t hat both the crite r ia by 
iMardling fe Aarsethl (120011) and IGiunnone et all (120131 ) 
are rather conservative because they have f e j, /star > 0.93 
and / 2 P i < 0.7. Thus, the systems satisfying these cri¬ 
teria are expected to be long-term stable, but those sys¬ 
tems that do not satisfy this condition are not necessarily 
expected to be unstable. 

Finally, we observe that the empirical stability bound¬ 
ary by LEggleton fe Kiseleval (|1995l ) performs the best 
among the previously proposed criteria. From Table 1, 
we observe that / 2 P i — 0.92, / e j ~ 0.84, and / star — 0.87, 
which are comparable to those obtained from our one- 
parameter criterion r ap = 1.83, but significantly lower 
than our two-parameter boundaries in Equations m 
and (fl6l) 

In summary, the_ stability boundary by 

lEggleton fc Kiseleval (ll995l) performs the best among 
the previously proposed crit eria, while those by 
IMardling fe Aarsethl ([20011 4 and IGiunnone et al.l ([2013( 4 
are too conservative. The Hill stability criterion has poor 
performance and provides very little useful information 
regarding the fate of the Hill unstable systems. 

5.2. Relation to other works with more than two planets 

Our results are strictly valid only for two-planet sys¬ 
tems. However, some of our main findings can still pro¬ 
vide useful information regarding the long-term stability 
in systems with more than two planets. 

First, we have found that the stability boundary de¬ 
pends on the eccentricities only through r ap = a ou t(l — 
e 0 ut)/ain(l+ein), which means that the relevant quantity 
to describe the stability is the distance between the peri- 
center of the outer planet and the apocenter of the inner 
planet. Moreover, we show that the relative orientation 
of the ellipses plays only a minor role in separating a sta¬ 
ble from unstable systems (see 44.1.41) . Consisten t with 
our results, the recent experiments bv iPu fe Wul (120151 ) 
show a similar dependence on the stability boundary in 
systems with seven planets. In their study the relevant 
quantity is the distance between the pericenter of the 

3 Consistent with the definition of Hill stability, we have checked 
that in all the unstable systems with r ap > V/t 1 ) 1 it is the outer 
(inner) planet the one that is ejected (collides with the star). Oth¬ 
erwise the planets would have had orbit crossing events. 
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outer planet and the apocenter of the inner planet for all 
the adjacent planets. 

Second, we find that the stability boundary changes 
~ 3 times more slowly with the maximum integration 
time for t max /P in y < 10 6 than in the range t max /P in j = 
10 6 —10 8 (see Figured]). We ob serve a similar behavior i n 
the numerical experiments bv iSmith fe Lissauerl (l2009tl . 
wi th more than two plan ets (see Figure 1 therein) and 
bv lChatteriee et all (120081) with three planets (see Figure 
29 therein). There, the slope of the separation between 
adjacent planets required for stability as a function of 
tmax/Pin,i decreases significantly after ~ 10 6 orbits of 
the innermost planet in the system. 

Previous studies generally parameterized the spacing 
required for stability in units of the mutual Hill radii 
Ph as K = (a nil t, ~ Qin)/PH = a + Mog(t max ) (e.g., 
iChambers et al.l 119961 : ISmith fc Lissauerl I2009H . which 
makes it hard to make a direct quantitative compari¬ 
son with our results in Equation (1151) . However, we can 
approximate our stability boundary in Equation (fTSl) for 
the limit of small eccentricities (or a ou t — ciin "C a ou t) 
and equal-mass planets (/q n = /i ou t) to write K oc 
b[ 3/(2/q n )] 1/3 log(f max ), where b is the coefficient we have 
obtained from our simulations and the Support Vector 
Machine algorithm. From our fits in Figure [H we find 
b = 0.021 for t max /P in ,i = 10 4 - 10 6 and b = 0.069 


for tmax/Pin,i = 10 6 — 10 8 . Since our simulations have 
an average mass ratios /!;„ = /lout = 10 , we de¬ 

rive K oc 2.4 log(f max ) for t max /Pi n ,i = 10 4 - 10 6 and 
K oc 0.79log(t max ) for f max /P in ,i = 10 6 - 10 s . 

We notice that the slope of 0.79 that we found for 
tmax/Pin,i = 10 6 — 10 8 falls in the ran ge of b ~ 0.7 — 1.3 
that w a s found in previou s stud ies bv ISmith fc Lissaued 
(|2009ll . I Punk et al.l d20lch . and IPu fc Wd (12015)1) . Our 
results also show that it is not possible to fit the bound¬ 
ary with a single linear fit in log (f max) since the slope b 
decreases as a function of time. This expectation is con¬ 
sistent with previous studies that predict l ower values 
of b as the integration time increases f max : iFunk et al.1 
(1201011 predict b ~ 1.3 for f m ax/Pin,i = 10 4 — 10 7 , 
ISmith fc Lissauerl (|2009D predict b ~ 1 for f m ax/Pin,i < 
10 s , and IPu fc Wul (I2015T) predict b ~ 0.7 for i m ax/Pin,i = 
10 7 - 10 9 

In summary, our results show that the stability bound¬ 
ary depends on the eccentricities only through the dis¬ 
tance between the orbits and logarithmically on the max¬ 
imum integration time, with a shallower slope for longer 
times. These results are qualitatively consistent with 
other stability studies with more than two planets. 


5.3. Application to observed planetary systems 

The results from our fiducial simulation 2pl-fiducial 
show that (see panel (b) in Figures 0] and [6j): 

1. with probability > 0.95 a system is unstable if 

I/O 

r ap < 2.4 [max(/ii n ,/tout)] ; +0.6, with ejections 
occurring for /q n > /i out and either ejection or col¬ 
lisions with the star for /q n < /x out ; 

2. and with probability > 0.95 a system is stable 
against either ejections or collisions with the star if 

r ap > 2.4 [max(/i in ,/tout )] 1/3 (aout/oin ) 1 - 72 + 1-4. 
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Fig. 9.— Stability boundary in Equation GH> as a function 
of the semi-major axis ratio a 0 ut/ain for a sample of two-planet 
syste ms discovered by rad ial velocity surveys wi th aout/am < 5 
(from I Wright et al.l 1201 U and HD 67087 from fHarakaw a et al.l 
2015). The error bars only consider the errors in the eccentricities 
and we use the minimum planet masses to calculate /i- in and // 0 ut- 
The vertical red dashed lines indicate the regions for which > 95% 
of the systems to the left (right) are expected to be unstable 
(stable) according to the stability criterion. The horizontal dashed 
lines indicate the position of the strongest mean-motion resonances. 


In Figure [9] we show the stability boundary in Equa¬ 
tion E] for a sample of two-planet systems discovered by 
radial velocity surveys. We also display the regions for 
which > 95% of the systems to the left (right) are ex¬ 
pected to be unstable (stable) according to the criterion 
above. 

We observe that some systems are expected to be un¬ 
stable according to our results. In particular, there are 3 
and 4 systems around the 2 : 1 and 3 : 2 that are consis¬ 
tent with being left to the dashed vertical, respectively. 
These results might seem to contradict the validity of our 
stability constraints. However, our results apply to more 
widely-spaced systems with a ou t/<Jin > 3, where we avoid 
the effect from these first-order mean-motion resonances 
that can promote the long-term stability of the system. 

A more curious case is the two-planet system HD 
202206 because it has a ou t/ain = 3.1 and our results 
should apply t o this rang e of a m it/a\n- A s discu ssed by 
iCorreia et al.l (|2005f ) and iCouetdic et al.l (|2010H such a 
system is indeed unstable for the best three-body fit of 
the RV measurements. However, there are stable copla- 
nar solutions provided that the system is in a 5 : 1 mean- 
motion reso nance. _ 

Recently, lHarakawa et al.l (12015T ) discovered the the 
planetary system HD 67087, which contains two planets 
with minimum masses of fi ln ~ 0.002 and p, out ~ 0.004 
and orbital elements a out /aj n ~ 3.6+ 9 ' 24 > e in = 0.171 qq 7 , 
and e out = OTO^q'^- This system is particularly inter¬ 
esting because our stability criterion indicates that the 
systems should be unstable unless the value of e ou t is in 
the lower end of its error measurement (see the error bar 
in Figure 9). This result suggests that the dynamical 
stability of this system should be further investigated, 
including the possibility of non-coplanar configurations 
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of the orbits that can lead to more stable solutions. 

In conclusion, all of the observed systems (with the 
exception of HD 67087) that are likely to be unstable 
according to our criterion seem to be protected by a 
mean-motion resonance. The effect of mean-motion res¬ 
onances does not play a significant role in our calcula¬ 
tions because we have excluded the lowest-order mean- 
motion resonances p : p + q with p > 0 and q = {1,2,3} 
(aout/ain = (1-58,2.08,2.51}) from our calculations. 

Finally, our results can be used to put constraints on 
the orbital elements of potential planets in systems with 
RV trends or poorly constrained RV measurements. In 
what follows, we give one worked example where we ap¬ 
ply our stability boundary. 

5.3.1. A worked example: constraints on the eccentricity of 

KOI-1299c. 

KOI-1299 is a giant star harboring at least two gi¬ 
ant planets (e.g., ICiceri et all 120151 IQrtiz et al.1 120151 
IQuinn et al.1 l2015h . The planetary system is in a hier¬ 
archical configuration with a ou t/ain — 4 and the inner 
planet (KOI-1299b) is in an eccentric orbit ei n ~ 0.5. The 
planet-to-mass ratios are p m ~ 0.004 and p ou t > 0.0018. 

Using these parameters and assuming that the plan¬ 
ets have relatively low mutual inclinations (i m < 40°) 
and p,[ n > /Tout, our stability constraint above implies 
that the system is unstable against ejections with prob¬ 
ability > 95% if the outer planet has an eccentricity of 
eout ^ 0.5. Therefore, we conclude that with high proba¬ 
bility that the eccentricity of KOI-1299c is e c < 0.5. This 
upper limit is useful in this example because the RV mea¬ 
surements by IQuinn et al.1 (1201511 yield e c = 0.64+°; 
and the error bar can be shrunk by using our stability 
constraint. Consistently, the authors have indeed studied 
the stability of this system and concluded that the system 
can be stable in a coplanar configuration for ~ 6 x 10 6 
orbits of the inner planet only if e c < 0.55, which then al¬ 
lowed them to fit the orbital parameters to much higher 
accuracy, finding e c = 0.498+ 0 ; 059 . 

It might be surp rising that the upper limit found by 
IQuinn et al. (120151 ) is less constraining than the one we 
found with our stability boundary. However, these au¬ 
thors study the stability of the system surveying a much 
more constrained region of parameter space, confining 
the orbits to be almost apsidally aligned, which allows 
for stable orbits with higher values of e c compared to 
random orientation of the orbits as we have assumed in 
our simulations (see H4.1.4I) . 

We repeated the analysis above by using all the 
others stability boundaries shown in Figure 0] to de¬ 
termine an upper limit to e c . The single-parameter 
bou ndary re qu ires e r S. 0. 59, while the boundary from 
lEggleton fc Kiseleval ( 1995 ) in Equation (j4j e c < 0.57. 
All other stability boundaries that we have tested here 
(panels (d), (e), and (f) in FigureQ} do not provide a use¬ 
ful constraint as they only demand e c < 1 for ejections 
not occur with probability > 0.95. 

In summary, our stability constraint places a strong 
constraint on the eccentricity of KOI -1299c, which iscon- 
sistent with the stability analysis of IQuinn et al. ( I2015h 
for this system. All other prev iousl y proposed s tability 
boundaries, except that of lEggleton fc Kiseleval (|1995l l. 
do not place a useful constraint to the eccentricity of 
KOI-1299c. 


5.4. Stellar evolution and white dwarf pollution 

From Table 1, we note that the ratio between the num¬ 
ber of stellar collisions and ejections is ~ 0.36, mean¬ 
ing that ~ 27% of the unstable systems reach distances 
< Rq if the inner planet starts a; n ,i = 46.5 AU. This 
fraction increases up to ~ 41% by placing the planet at 
din,i = 0.465 AU (see Figure 0 ■ Since a Jupiter-like 
planet orbiting a 0 . 6 M 0 white dwarf is expected to be 
disrupted in a highly eccentric orbit i f it reaches a dis¬ 
tance < 3 Rq (iGuillochon et al.l I 2011 H . we expect that 
unstable hierarchical two-planet systems can often lead 
to tidal disruptions. 

Note that most (~ 95%) stellar collisions start with 
Sin < Sout (see Figure 0. Also, we observe that the 
ratio between the number stellar collisions and the num¬ 
ber of ejections in our fiducial simulation is < 0.1 for 
Mout/p-in 1-5 and it reaches values of ~ 1 — 3 for 
/Xout/ltin ~ 2 — 6 . These results are qualitatively consis¬ 
tent with the increase in the ratio between the number 
of planets undergoing a close approach with the star and 
the number of ejections from ~ 0.03 for equal-mass plan¬ 
ets to ~ 0.12 —0.16 for planetary-mass ratios of ~ 2.3 — 3 
(randomly assigning the more mass ive planet as the in¬ 
ner one) observed bv lFord fc Rasiol (l2008ll . However, we 
observe that the overall rate of collisions with the star 
relative to ejections can be several times higher in our 
simulations for two initiall y eccentric p lanets relative to 
the simulations by iFord fc Rasiol (12008H for two planets 
in initially circular orbits. 

Equation (fl5|) shows that as the planetary system ages 
the our stability boundary becomes more stringent, al¬ 
lowing orbits with relatively larger separations (larger 
r ap ) to become unstable. However, the dependence is 
only logarithmic and the boundary moves only by ~ 7% 
percent per order magnitude difference in the evolution 
time. Thus, by extrapolating this result to timescales 
> 10 8 Pi n one would expect only a small effect in the 
stability of planetary systems. 

A more pronounced effect from the aging of the plan¬ 
etary system is likel y to co me from mass loss of the 
host star fe.g. JDebes fe Sigurdssonir2002D . Typical white 
dwarfs have masses that are a few times lower than their 
main-sequence progenitors and therefore the mass ratios 
Si n and p out are expected to increase by th e same fac¬ 
tor, while keeping a n ut/a\r, fixed (see, iVeras et al.l 120131 
iMustill et al.l 12014 : iVeras fc Gansickdl2015fl . This effect 
is expected to destabilize the systems close to our stabil¬ 
ity boundary in Equation m- 

In summary, unstable two-planet systems in an ini¬ 
tially hierarchical configuration can lead to a significant 
number of collisions with the star relative to the number 
of ejections, which might contribute to the pollution of 
white dwarfs as a result of stellar mass loss. The number 
of collisions with the star (or tidal disruptions) can be 
higher than the number of ejections for /lout/Min ~ 2 — 6 . 

6 . CONCLUSIONS 

We run a large number of long-term numerical inte¬ 
grations to study the fates of two-planet systems in hier¬ 
archical configurations with arbitrary eccentricities and 
mutual inclinations. 

Using the Support Vector Machine algorithm to sep¬ 
arate different fates of our simulated systems, we 
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find that initially nearly coplanar systems remain 
long-term stable for a ou t(l - e ou t)/[ain(l + e in )] > 

2.4 [max(/i in , /x out )] 1/3 (a ou t/ai„) 1/2 + 1-15. Systems that 
do not satisfy this condition by a margin of > 0.5 are ex¬ 
pected to be unstable, mostly leading to planet ejections 
if /iin > /Tout, while slightly favoring collisions with the 
star for /x in < /x out . 

We show that our proposed stability boundary 
performs significantly bet ter th an previously pro¬ 
posed stab ility _criteria (lEggleton &; Kiseleva! I1995L 
iMardling fe Aarsethll200il and Hill stability) for mutual 
inclinations < 40°. 
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